rm(list=ls())
library(foreign)
library(psy) #Calcuate Cronbach's alphas.
library(RI2by2) #Randomization inference for binary outcome variables.
library(ipumsr) #Import population values from the IPUMS data.
library(tidyverse)
library(dplyr)


####Import data.

####Sample
###Sample Size
###Cronbach's Alpha for Composite Indices 
###Table01_Balance: Checks for Balance
####Results
###Table02_Democracy: Table for the Effect of Democracy
###Table03_DemocracyUN: Table for the Effect of Democracy, Controlling for UN Approval
####Robustness Checks
###Table04_Logit: Table of Logit Regressions with demo*un Interaction 
###Checks for Heterogeneity Treatment Effects
###TableA09_Specific: Breakdown of the Open-Ended Responses of Non-Compliers
###Figure 1: Predicted Probabilities from Logistic Regressions
###Plug In Estimator
###Logistic Regressions and Predicted Probabilities with Population Values
###TableA10_DemocracyControls: Effect of Democracy, Controlling for Attitudinal and Demographic Attiributes
###Randomization Inference


####Read the data files into R
setwd("~/Dropbox/AudCost/RnR")
Brazil<-read.csv("Brazil.csv", header=TRUE)
China<-read.csv("China.csv", header=TRUE)
YouGov<-read.csv("YouGov.csv", header=TRUE)

##Dataframes With Only Complete Entries for Predicted Probabilities and Post-Stratification
myvars <- c("cow_ccode", "attack_support", "vignette", "democracy", "un", "militarism", "internationalism", "nationalism", "specific", "age", "female", "education", "education2", "income", "news", "religious", "religiosity")
Brazil_n <- na.omit(Brazil[myvars])
China_n <- na.omit(China[myvars])
YouGov_n<-na.omit(YouGov[myvars])


####Sample

###Sample Size
nrow(Brazil)
nrow(China)
nrow(YouGov)

###Cronbach's Alpha for Composite Indices 
Militarism_Brazil<-cbind (Brazil$militarism_1, Brazil$militarism_2)
cronbach(Militarism_Brazil)

Internationalism_Brazil<-cbind (Brazil$internationalism_1, Brazil$internationalism_2, Brazil$internationalism_3, Brazil$internationalism_4)
cronbach(Internationalism_Brazil)

Nationalism_Brazil<-cbind (Brazil$nationalism_1, Brazil$nationalism_2)
cronbach(Nationalism_Brazil)

Militarism_China<-cbind (China$militarism_1, China$militarism_2)
cronbach(Militarism_China)

Internationalism_China<-cbind (China$internationalism_1, China$internationalism_2, China$internationalism_3, China$internationalism_4)
cronbach(Internationalism_China)

Nationalism_China<-cbind (China$nationalism_1, China$nationalism_2)
cronbach(Nationalism_China)

###Table01_Balance: Checking for Balance

##Descriptive Statistics
T1<-function(FUN,d,v){
  m<-FUN(d[[v]], na.rm=TRUE)
  return(m)
}

T1(mean, Brazil, "militarism")
T1(mean, China, "militarism")
T1(mean, YouGov, "militarism")

T1(mean, Brazil, "internationalism")
T1(mean, China, "internationalism")
T1(mean, YouGov, "internationalism")

T1(mean, Brazil, "nationalism")
T1(mean, China, "nationalism")
T1(mean, YouGov, "nationalism")

T1(mean, Brazil, "religiosity")
T1(mean, China, "religiosity")
T1(mean, YouGov, "religiosity")

T1(mean, Brazil, "age")
T1(mean, China, "age")
T1(mean, YouGov, "age")

T1(mean, Brazil, "news")
T1(mean, China, "news")
T1(mean, YouGov, "news")

T1(median, Brazil, "education")
T1(median, China, "education")
T1(median, YouGov, "education")

T1(median, Brazil, "income")
T1(median, China, "income")
T1(median, YouGov, "income")

prop.table(table(Brazil$female))
prop.table(table(China$female))
prop.table(table(YouGov$female))

prop.table(table(Brazil$religious))
prop.table(table(China$religious))
prop.table(table(YouGov$religious))

nrow(Brazil)
nrow(China)
nrow(YouGov)

###Checking for Balance for the Brazil Sample
anova(lm(age ~ vignette, Brazil))
anova(lm(education ~ vignette, Brazil))
anova(lm(income ~ vignette, Brazil))
anova(lm(female ~ vignette, Brazil))
anova(lm(news ~ vignette, Brazil))
anova(lm(religious ~ vignette, Brazil)) 
anova(lm(religiosity ~ vignette, Brazil)) #Unbalanced at the 0.05 level

###Checking for Balance for the China Sample
anova(lm(age ~ vignette, China))
anova(lm(education ~ vignette, China))
anova(lm(income ~ vignette, China)) 
anova(lm(female ~ vignette, China))
anova(lm(news ~ vignette, China)) 
anova(lm(religious ~ vignette, China)) 
anova(lm(religiosity ~ vignette, China))

###Checking for Balance for the YouGov Sample
anova(lm(age ~ vignette, YouGov))
anova(lm(education ~ vignette, YouGov))
anova(lm(income~ vignette, YouGov)) 
anova(lm(female ~ vignette, YouGov))
anova(lm(news ~ vignette, YouGov)) 
anova(lm(religious ~ vignette, YouGov)) 
anova(lm(religiosity ~ vignette, YouGov))


####Results

###Table02_Democracy: Effect of Democracy
tapply(Brazil$attack_support, list("Democracy"=Brazil$democracy), mean, na.rm=TRUE)
summary(lm(attack_support ~ democracy, data = Brazil))

tapply(China$attack_support, list("Democracy"=China$democracy), mean, na.rm=TRUE)
summary(lm(attack_support ~ democracy, data = China)) 

tapply(YouGov$attack_support, list("Democracy"=YouGov$democracy), mean, na.rm=TRUE)
summary(lm(attack_support ~ democracy, data = YouGov)) 


###Table03_DemocracyUN: Effect of Democracy, Controlling for UN Approval

##Effect of Democracy
tapply(Brazil$attack_support, list("Democracy"=Brazil$democracy, "UN Approval"=Brazil$un), mean, na.rm=TRUE)
summary(lm(attack_support ~ democracy, data = Brazil[Brazil$un==1,])) 
summary(lm(attack_support ~ democracy, data = Brazil[Brazil$un==0,])) 

tapply(China$attack_support, list("Democracy"=China$democracy, "UN Approval"=China$un), mean, na.rm=TRUE)
summary(lm(attack_support ~ democracy, data = China[China$un==1,])) 
summary(lm(attack_support ~ democracy, data = China[China$un==0,])) 

tapply(YouGov$attack_support, list("Democracy"=YouGov$democracy, "UN Approval"=YouGov$un), mean, na.rm=TRUE)
summary(lm(attack_support ~ democracy, data = YouGov[YouGov$un==1,])) 
summary(lm(attack_support ~ democracy, data = YouGov[YouGov$un==0,])) 

##Interaction democracy*un is not significant.
summary(lm(attack_support ~ democracy*un, data = Brazil)) 
summary(lm(attack_support ~ democracy*un, data = China)) 
summary(lm(attack_support ~ democracy*un, data = YouGov)) 
                        
##Effect of UN Approval
summary(lm(attack_support ~ un, data = Brazil[Brazil$democracy==1,])) 
summary(lm(attack_support ~ un, data = Brazil[Brazil$democracy==0,])) 

summary(lm(attack_support ~ un, data = China[China$democracy==1,]))
summary(lm(attack_support ~ un, data = China[China$democracy==0,]))

summary(lm(attack_support ~ un, data = YouGov[YouGov$democracy==1,])) 
summary(lm(attack_support ~ un, data = YouGov[YouGov$democracy==0,])) 


###Table04_Logit: Logit Regressions with Democracy*UN Interaction 

##Logistic Regressions
BrazilLogitGLM<-glm(attack_support ~ democracy*un + militarism + internationalism + nationalism + specific
                    + age + female + education + income + news + religious + religiosity, 
                    family = "binomial", data = Brazil_n)
summary(BrazilLogitGLM)

ChinaLogitGLM<-glm(attack_support ~ democracy*un + militarism + internationalism + nationalism + specific
                   + age + female + education + income + news + religious + religiosity, 
                   family = "binomial", data = China_n)
summary(ChinaLogitGLM)

YouGovLogitGLM<-glm(attack_support ~ democracy*un + militarism + internationalism + nationalism + specific
                    + age + female + education + income + news + religious + religiosity, 
                    family = "binomial", data = YouGov_n)
summary(YouGovLogitGLM)

###Difference Among Samples
mean(Brazil$age)
sd(Brazil$age)

mean(China$age)
sd(China$age)

mean(YouGov$age)
sd(YouGov$age)

100*prop.table(table(Brazil$education))
100*prop.table(table(China$education))

100*prop.table(table(Brazil$religious))
100*prop.table(table(China$religious))

mean(Brazil$militarism, na.rm=TRUE)<mean(China$militarism, na.rm=TRUE)

mean(Brazil$nationalism, na.rm=TRUE)<mean(China$nationalism, na.rm=TRUE)

### Heterogeneity by context?
## Interaction of thinking of specific case and democracy not significant
summary(lm(attack_support ~ democracy*specific, data = Brazil))

## Which country you ID does not affect ATE
m1<-lm(attack_support ~ democracy*open_end, data = Brazil)
m2<-lm(attack_support ~ democracy+open_end, data = Brazil)
anova(m1,m2)

## Interaction of thinking of specific case and democracy not significant
summary(lm(attack_support ~ democracy*specific, data = China))

## Which country you ID does not affect ATE
m3<-lm(attack_support ~ democracy*open_end, data = China)
m4<-lm(attack_support ~ democracy+open_end, data = China)
anova(m3,m4)


###TableA09_Specific: breakdown of the open-ended responses of non-compliers
sort(table(Brazil$open_end), decreasing=TRUE)[1:11]
sort(100*table(Brazil$open_end)/length(Brazil$specific[Brazil$specific==1]), decreasing=TRUE)[1:11]
length(Brazil$specific[Brazil$specific==1])
length(Brazil$specific[Brazil$specific==0])

sort(table(China$open_end), decreasing=TRUE)[1:11]
sort(100*table(China$open_end)/length(China$specific[China$specific==1]), decreasing=TRUE)[1:11]
length(China$specific[China$specific==1])
length(China$specific[China$specific==0])

sort(table(YouGov$open_end), decreasing=TRUE)[1:11]
sort(100*table(YouGov$open_end)/length(YouGov$specific[YouGov$specific==1]), decreasing=TRUE)[1:11]
length(YouGov$specific[YouGov$specific==1])
length(YouGov$specific[YouGov$specific==0])


###Figure 1: Predicted Probabilities from Logistic Regressions
T4<-function(d, x, y, m){
  vig <- with(d, data.frame(democracy=x,
                             un=y,
                             militarism=mean(d[["militarism"]]), 
                             internationalism=mean(d[["internationalism"]]), 
                             nationalism=mean(d[["nationalism"]]), 
                             specific=median(d[["specific"]]), 
                             age=mean(d[["age"]]), 
                             female=median(d[["female"]]), 
                             education=median(d[["education"]]), 
                             income=median(d[["income"]]), 
                             news=mean(d[["news"]]), 
                             religious=median(d[["religious"]]), 
                             religiosity=mean(d[["religiosity"]])))
  preds <- data.frame(predict(m, vig, type="response", se.fit=TRUE))
  return(preds)
}

T4(Brazil_n, 0, 1, BrazilLogitGLM)
T4(Brazil_n, 1, 1, BrazilLogitGLM)
T4(Brazil_n, 0, 0, BrazilLogitGLM)
T4(Brazil_n, 1, 0, BrazilLogitGLM)

T4(China_n, 0, 1, ChinaLogitGLM)
T4(China_n, 1, 1, ChinaLogitGLM)
T4(China_n, 0, 0, ChinaLogitGLM)
T4(China_n, 1, 0, ChinaLogitGLM)

T4(YouGov_n, 0, 1, YouGovLogitGLM)
T4(YouGov_n, 1, 1, YouGovLogitGLM)
T4(YouGov_n, 0, 0, YouGovLogitGLM)
T4(YouGov_n, 1, 0, YouGovLogitGLM)

###Plug In Estimator
democracy_B2<-Brazil_n$democracy+1
by(BrazilLogitGLM$fitted.values,democracy_B2,mean)
by(BrazilLogitGLM$fitted.values,list(democracy_B2,Brazil_n$un),mean)

democracy_C2<-China_n$democracy+1
by(ChinaLogitGLM$fitted.values,democracy_C2,mean)
by(ChinaLogitGLM$fitted.values,list(democracy_C2,China_n$un),mean)

democracy_Y2<-YouGov_n$democracy+1
by(YouGovLogitGLM$fitted.values,democracy_Y2,mean)
by(YouGovLogitGLM$fitted.values,list(democracy_Y2,YouGov_n$un),mean)


###Logistic Regressions and Predicted Probabilites with Population Values

##Logistic Regressions
BrazilLogitGLM_2<-glm(attack_support ~ democracy*un + militarism + internationalism + nationalism + specific
                    + age + female + education2 + income + news + religious + religiosity, 
                    family = "binomial", data = Brazil_n)

ChinaLogitGLM_2<-glm(attack_support ~ democracy*un + militarism + internationalism + nationalism + specific
                   + age + female + education2 + income + news + religious + religiosity, 
                   family = "binomial", data = China_n)

YouGovLogitGLM_2<-glm(attack_support ~ democracy*un + militarism + internationalism + nationalism + specific
                     + age + female + education2 + income + news + religious + religiosity, 
                     family = "binomial", data = YouGov_n)

##Download individual-level data from the IPUMS website (https://international.ipums.org/international/)
##In particular, we used the variables "AGE", "AGE2", "SEX", "EDATTAIN" from Brazil's 2010 Census and China's 2000 Census.
##Note: we do not have the authority to distribute the data to a third party; the following R code was provided by the IPUMS. 
setwd("~/Dropbox/IPUMS_2018")
if (!require("ipumsr")) stop("Reading IPUMS data into R requires the ipumsr package. It can be installed using the following command: install.packages('ipumsr')")
ddi <- read_ipums_ddi("ipumsi_00008.xml")
ipums_data <- read_ipums_micro(ddi)
ipums_data<-
  ipums_data %>%
  select(COUNTRY, AGE, AGE2, SEX, EDATTAIN) %>%
  mutate(EDATTAIN = na_if(EDATTAIN, 0), #Replace NIUs with NAs.
         female=ifelse(SEX==1, 0, 1), 
         education2=ifelse(EDATTAIN==1, 1, EDATTAIN-1)
         ) %>% 
  filter(AGE>=18) %>% #Drop minors.
  dplyr::rename(age=AGE) %>%
  select(COUNTRY, age, female, education2) %>%
  drop_na()

Brazil_pop<-ipums_data[ipums_data$COUNTRY==76,]
China_pop<-ipums_data[ipums_data$COUNTRY==156,]

T4_p<-function(d, p, x, y, m){
  vig <- with(d, data.frame(democracy=x,
                            un=y,
                            militarism=mean(d[["militarism"]]), 
                            internationalism=mean(d[["internationalism"]]), 
                            nationalism=mean(d[["nationalism"]]), 
                            specific=median(d[["specific"]]), 
                            age=mean(p[["age"]]), 
                            female=median(d[["female"]]), 
                            education2=median(p[["education2"]]), 
                            income=median(d[["income"]]), 
                            news=mean(d[["news"]]), 
                            religious=median(d[["religious"]]), 
                            religiosity=mean(d[["religiosity"]])))
  preds <- data.frame(predict(m, vig, type="response", se.fit=TRUE))
  return(preds)
}

T4_p(Brazil_n, Brazil_pop, 0, 1, BrazilLogitGLM_2)
T4_p(Brazil_n, Brazil_pop, 1, 1, BrazilLogitGLM_2)
T4_p(Brazil_n, Brazil_pop, 0, 0, BrazilLogitGLM_2)
T4_p(Brazil_n, Brazil_pop, 1, 0, BrazilLogitGLM_2)

T4_p(China_n, China_pop, 0, 1, ChinaLogitGLM_2)
T4_p(China_n, China_pop, 1, 1, ChinaLogitGLM_2)
T4_p(China_n, China_pop, 0, 0, ChinaLogitGLM_2)
T4_p(China_n, China_pop, 1, 0, ChinaLogitGLM_2)

T4_p(YouGov_n, China_pop, 0, 1, YouGovLogitGLM_2)
T4_p(YouGov_n, China_pop, 1, 1, YouGovLogitGLM_2)
T4_p(YouGov_n, China_pop, 0, 0, YouGovLogitGLM_2)
T4_p(YouGov_n, China_pop, 1, 0, YouGovLogitGLM_2)


###TableA10_DemocracyControls

TA10_1<-function(d,v,l){
  if(l=="L"){
  attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[[v]]==0]
  attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[[v]]==0]
  }
  if(l=="M"){
  attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[[v]]==0.5]
  attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[[v]]==0.5]
  }
  if(l=="H"){
  attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[[v]]==1]
  attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[[v]]==1]
  }
  return(t.test(attack_support.auto,attack_support.demo))
}

TA10_1(Brazil, "militarism", "L")
TA10_1(Brazil, "militarism", "M")
TA10_1(Brazil, "militarism", "H")

TA10_1(China, "militarism", "L")
TA10_1(China, "militarism", "M")
TA10_1(China, "militarism", "H")

TA10_intlism_natlism<-function(d,v,l){
  if(l=="L"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[[v]]<=0.25]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[[v]]<=0.25]
  }
  if(l=="M"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & 0.25<d[[v]] & d[[v]]<0.75]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & 0.25<d[[v]] & d[[v]]<0.75]
  }
  if(l=="H"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[[v]]>=0.75]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[[v]]>=0.75]
  }
  return(t.test(attack_support.auto,attack_support.demo))
}

TA10_intlism_natlism(Brazil, "internationalism", "L")
TA10_intlism_natlism(Brazil, "internationalism", "M")
TA10_intlism_natlism(Brazil, "internationalism", "H")

TA10_intlism_natlism(China, "internationalism", "L")
TA10_intlism_natlism(China, "internationalism", "M")
TA10_intlism_natlism(China, "internationalism", "H")

TA10_intlism_natlism(Brazil, "nationalism", "L")
TA10_intlism_natlism(Brazil, "nationalism", "M")
TA10_intlism_natlism(Brazil, "nationalism", "H")

TA10_intlism_natlism(China, "nationalism", "L")
TA10_intlism_natlism(China, "nationalism", "M")
TA10_intlism_natlism(China, "nationalism", "H")

TA10_1(Brazil, "specific", "L")
TA10_1(Brazil, "specific", "H")

TA10_1(China, "specific", "L")
TA10_1(China, "specific", "H")

TA10_1(Brazil, "female", "H")
TA10_1(Brazil, "female", "L")

TA10_1(China, "female", "H")
TA10_1(China, "female", "L")

TA10_1(Brazil, "religious", "H")
TA10_1(Brazil, "religious", "L")

TA10_1(China, "religious", "H")
TA10_1(China, "religious", "L")

TA10_education<-function(d,l){
  if(l=="less than high school diploma"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & (d[["education"]]==1|d[["education"]]==2)]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & (d[["education"]]==1|d[["education"]]==2)]
  }
  if(l=="high school diploma"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["education"]]==3]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["education"]]==3]
  }
  if(l=="some college (degree)"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & (d[["education"]]==4|d[["education"]]==5)]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & (d[["education"]]==4|d[["education"]]==5)]
  }
  if(l=="some graduate school (degree)"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & (d[["education"]]==6|d[["education"]]==7)]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & (d[["education"]]==6|d[["education"]]==7)]
  }
  return(t.test(attack_support.auto,attack_support.demo))
}

TA10_education(Brazil, "less than high school diploma")
TA10_education(Brazil, "high school diploma")
TA10_education(Brazil, "some college (degree)")
TA10_education(Brazil, "some graduate school (degree)")

TA10_education(China, "less than high school diploma")
TA10_education(China, "high school diploma")
TA10_education(China, "some college (degree)")
TA10_education(China, "some graduate school (degree)")

TA10_income<-function(d,l){
  if(l=="q12"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & (d[["income"]]==1|d[["income"]]==2)]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & (d[["income"]]==1|d[["income"]]==2)]
  }
  if(l=="q3"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["income"]]==3]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["income"]]==3]
  }
  if(l=="q4"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["income"]]==4]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["income"]]==4]
  }
  if(l=="q5"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["income"]]==5]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["income"]]==5]
  }
  return(t.test(attack_support.auto,attack_support.demo))
}

TA10_income(Brazil, "q12")
TA10_income(Brazil, "q3")
TA10_income(Brazil, "q4")
TA10_income(Brazil, "q5")

TA10_income(China, "q12")
TA10_income(China, "q3")
TA10_income(China, "q4")
TA10_income(China, "q5")

TA10_news<-function(d,l){
  if(l=="L"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["news"]]<=2]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["news"]]<=2]
  }
  if(l=="M"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & 2<d[["news"]] & d[["news"]]<6]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & 2<d[["news"]] & d[["news"]]<6]
  }
  if(l=="H"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["news"]]>=6]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["news"]]>=6]
  }
  return(t.test(attack_support.auto,attack_support.demo))
}

TA10_news(Brazil, "L")
TA10_news(Brazil, "M")
TA10_news(Brazil, "H")

TA10_news(China, "L")
TA10_news(China, "M")
TA10_news(China, "H")

TA10_religiosity<-function(d,l){
  if(l=="L"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["religiosity"]]==0]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["religiosity"]]==0]
  }
  if(l=="M"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & 0<d[["religiosity"]] & d[["religiosity"]]<0.75]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & 0<d[["religiosity"]] & d[["religiosity"]]<0.75]
  }
  if(l=="H"){
    attack_support.auto<-d[["attack_support"]][d[["democracy"]]==0 & d[["religiosity"]]>=0.75]
    attack_support.demo<-d[["attack_support"]][d[["democracy"]]==1 & d[["religiosity"]]>=0.75]
  }
  return(t.test(attack_support.auto,attack_support.demo))
}

TA10_religiosity(Brazil, "L")
TA10_religiosity(Brazil, "M")
TA10_religiosity(Brazil, "H")

TA10_religiosity(China, "L")
TA10_religiosity(China, "M")
TA10_religiosity(China, "H")


###Randomization Inference
##Randomization inference about the democracy treatment
mat_Brazil = matrix(c(table(Brazil$democracy,Brazil$attack_support)[2,2],
                     table(Brazil$democracy,Brazil$attack_support)[2,1],
                     table(Brazil$democracy,Brazil$attack_support)[1,2],
                     table(Brazil$democracy,Brazil$attack_support)[1,1]),2,2,byrow=TRUE)

set.seed(1234567) 
ae_CI_Brazil<-AE.CI(mat_Brazil,0.05)
ae_CI_Brazil

mat_China = matrix(c(table(China$democracy,China$attack_support)[2,2],
                    table(China$democracy,China$attack_support)[2,1],
                    table(China$democracy,China$attack_support)[1,2],
                    table(China$democracy,China$attack_support)[1,1]),2,2,byrow=TRUE)
set.seed(1234567) 
ae_CI_China<-AE.CI(mat_China,0.05)
ae_CI_China

mat_YouGov = matrix(c(table(YouGov$democracy,YouGov$attack_support)[2,2],
                     table(YouGov$democracy,YouGov$attack_support)[2,1],
                     table(YouGov$democracy,YouGov$attack_support)[1,2],
                     table(YouGov$democracy,YouGov$attack_support)[1,1]),2,2,byrow=TRUE)
set.seed(1234567) 
ae_CI_YouGov<-AE.CI(mat_YouGov,0.05)
ae_CI_YouGov

##Randomization inference about the UN treatment
mat_Brazil2 = matrix(c(table(Brazil$un,Brazil$attack_support)[2,2],
                      table(Brazil$un,Brazil$attack_support)[2,1],
                      table(Brazil$un,Brazil$attack_support)[1,2],
                      table(Brazil$un,Brazil$attack_support)[1,1]),2,2,byrow=TRUE)

set.seed(1234567) 
ae_CI_Brazil2<-AE.CI(mat_Brazil2,0.05)
ae_CI_Brazil2

mat_China2 = matrix(c(table(China$un,China$attack_support)[2,2],
                     table(China$un,China$attack_support)[2,1],
                     table(China$un,China$attack_support)[1,2],
                     table(China$un,China$attack_support)[1,1]),2,2,byrow=TRUE)
set.seed(1234567) 
ae_CI_China2<-AE.CI(mat_China2,0.05)
ae_CI_China2

mat_YouGov2 = matrix(c(table(YouGov$un,YouGov$attack_support)[2,2],
                      table(YouGov$un,YouGov$attack_support)[2,1],
                      table(YouGov$un,YouGov$attack_support)[1,2],
                      table(YouGov$un,YouGov$attack_support)[1,1]),2,2,byrow=TRUE)
set.seed(1234567) 
ae_CI_YouGov2<-AE.CI(mat_YouGov2,0.05)
ae_CI_YouGov2

####THE END